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Analysis of Modulated Multivariate Oscillations 

Jonathan M. Lilly, Member, IEEE, and Sofia C. Olhede, Member, IEEE 



Abstract — The concept of a common modulated oscillation 
spanning multiple time series is formalized, a method for the 
recovery of such a signal from potentially noisy observations is 
proposed, and the time-varying bias properties of the recovery 
method are derived. The method, an extension of wavelet ridge 
analysis to the multivariate case, identifies the common oscillation 
by seeking, at each point in time, a frequency for which a 
bandpassed version of the signal obtains a local maximum in 
power. The lowest-order bias is shown to involve a quantity, 
termed the instantaneous curvature, which measures the strength 
of local quadratic modulation of the signal after demodulation 
by the common oscillation frequency. The bias can be made to 
be small if the analysis filter, or wavelet, can be chosen such that 
the signal's instantaneous curvature changes little over the filter 
time scale. An application is presented to the detection of vortex 
motions in a set of freely-drifting oceanographic instruments 
tracking the ocean currents. 

Index Terms — Amplitude and Frequency Modulated Signal, 
Analytic Signal, Bedrosian's Theorem, Complex- Valued Signal, 
Complex- Valued Time Series, Multivariate Signal. 

I. Introduction 

IN the physical sciences the description of common vari- 
ability in a set of multiple time series is an important 
data analysis task. Frequently a key signal present in the 
data is that of a modulated oscillation, extending across time 
series channels, with a different amplitude and perhaps a 
different phase shift in each time series. Such oscillations 
may be the signature of waves or wavelike phenomenon. 
The most important multivariate cases are the bivariate and 
trivariate cases, which occur frequently in oceanography and 
seismology, for example. One wishes to extract the common 
oscillatory structure from the observations, a task that is 
complicated by the possible presence of noise and also by 
time variability of the signal of interest. 

The analysis of univariate modulated oscillations is more 
highly evolved than is the multivariate case. In both cases 
one must begin with a model for the signal structure. For 
univariate signals, an attractive representation for an ampli- 
tude/frequency modulated signal is now well known, and 
involves the construction of a complex-valued quantity called 
the analytic signal 0]]-||6]]. Real-world signals are nearly 
always contaminated by noise or other sources of variability, 
hence some means of filtering or localizing the time series 
is required in order to isolate the modulated oscillation. The 
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analytic signal corresponding to the modulated oscillation of 
interest can be estimated with a popular and powerful method 
known as wavelet ridge analysis Q-O. The essence of this 
method, which is more general than its name might suggest, 
is a local optimization applied to a set of frequency-localized 
versions of the observed time series. 

This paper develops a powerful and flexible method, termed 
multivariate wavelet ridge analysis, for the extraction of 
modulated oscillations from multivariate time series. Estimates 
of the time-varying forms of leading-order bias terms are 
derived, which are essential in informing the choice of analysis 
filter or wavelet. This a non-trivial extension of a related work 
by the authors for the univariate case [0. The key innovation 
here is a model for signal structure in which a set of signals are 
expanded in terms of deviations from oscillatory behavior at a 
single common but time-varying frequency. The basic idea of 
multivariate wavelet ridge analysis, but without a theoretical 
understanding of the bias, was presented in the preliminary 
work ifTOl . 

The motivation for such a method is the analysis of ocean 
currents in the now very large set of data from freely- 
drifting, or "Lagrangian", instruments, see e.g. ifPD and ref- 
erences therein. The signatures of a particular type of oceanic 
structure — long-lived or "coherent" vortices ifTZl — occur fre- 
quently in such data and are aptly described as modulated 
oscillations in two dimensions. The development of automated 
and objective schemes for the analysis of such features has 
been attempted by several authors lfT3l - lfT5l . and thus this 
work will be of practical value. An application to a dataset of 
this type, from the observational experiment of |[T6l , IfPTl , is 
presented here as an illustration. 

The structure of the paper is as follows. Some essential 
background is presented in Section [II] together with a data 
example. In Section [Til] we introduce a representation for a 
modulated multivariate oscillation, and quantify the degree of 
variability of such a signal via a local expansion. A general- 
ization of wavelet ridge analysis appropriate to a multivariate 
signal is presented in Section [IV] and the leading-order bias 
term is identified. A key contribution is the identification and 
interpretation of the quantity controlling the bias, a higher- 
order relative of the joint instantaneous bandwidth of [18] 
which we term the joint instantaneous curvature. 

All data, numerical algorithms, and functions for analysis 
and figure generation are distributed to the community as a 
freely available Matlab package, as described in Appendix [aQ 

II. Background 

This section presents the background necessary for the 
development of an analysis method for treating modulated 

'This package, called Jlab, is available at http://www.jmlilly.net 
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multivariate oscillations. A real-world data example of oceano- 
graphic data provides a practical motivation. 



A. Statement of the Problem 

A set of N real-valued observed time series, assumed 
square-integrable herein, are arranged as an TV-vector 



Xo(t) = [x -i{t) x . 2 (t) ... x 0]N (t)] J 



(1) 



where "T" denotes the matrix transpose. At least some of 
the N channels of x D (f) are expected to contain oscillatory 
variability, and these oscillations are in turn expected to be 
related to one another or to share some joint structure. We 
therefore model x D (i) as containing two separate components 



x (i) = x(i) +x r (t) 



(2) 



where x(t) is a modulated multivariate oscillation, defined 
subsequently, and x r (t) is a residual which we assume may 
be accurately represented as a stochastic process. Thus x(t) 
is the "signal" and x r (t) is the "noise". Our goals here are 
(i) to estimate the multivariate oscillatory signal x(t) given 
the observed vector x Q (t); (ii) to characterize its time-varying 
behavior; and to estimate the errors in this process from (iii) 
bias associated with x(t) itself. The first step is a model 
specification for the modulated multivariate oscillation. 



B. A Bivariate Example 

An example of data matching the model (O for the bivariate 
case of N — 2 is shown in Fig. Q] along with our eventual 
decomposition into an estimated oscillatory portion x(t) plus 
an estimated residual x r (i). This data, described in more detail 
in Section [V] is from a set of freely-drifting instruments called 
"floats" that track the ocean currents, recording their horizontal 
position at regular intervals. Freely-drifting instruments such 
as these represent one of the primary ways oceanographers 
study the structure and variability of ocean currents, see e.g. 
ifTTI and references therein. 

The observed time series x a (t) in Fig. QJ clearly show 
the presence of modulated oscillations superposed on a 
background of apparently random fluctuations, matching the 
model (0. The oscillations in these records represent the 
presence of long-lived, intense oceanic vortex structures, in 
this case of about 10-20 km radius |fl6l , IfTTI . Vortices such as 
those seen here are ubiquitous features of the ocean currents 
[e.g., 12], and a large number of papers have been devoted 
to the study of their dynamics and impact on the large-scale 
flow. In Fig. [Th, the estimated bivariate modulated oscillations 
x(t) have been visually represented as time-varying ellipses; 
we refer the reader to lfj"8l for details on the modulated 
elliptical signal representation of a bivariate analytic signal. 
Such modulated elliptical signals are a special case of a more 
general class of signals we will consider; it is worth pointing 
out that ellipses form the building blocks for models of many 
different kinds of data, from ocean currents |fl9l to seismic 
signals ll20l to electroencephalographic (EEG) data lETI . 



C. Fundamentals 

The starting point for our analysis is the analytic signal 
method [ 1 1-| 6 1 for assigning meaningful time-varying ampli- 
tudes and frequencies to each of the channels of x(t). The 
analytic part of the signal vector is defined as 

x+(f) = 24x(t) = x(f) + iHx(t) (3) 
where '"H" denotes the Hilbert transform operator 

x(r) 



1 f 

Hxlt) = -4- 



t- 



dr 



(4) 



with being the Cauchy principal value integral. The 
Fourier transforms of x(i) and x + (t) are denoted X(w) and 
X+(w), respectively. It follows from the form of the analytic 
operator 4 in the frequency domain that 

X+(w) = 2U(u)X.(w) (5) 

where U (lo) is the unit step function. Thus the application of 
the operator 24 to x(t) doubles the amplitudes of the Fourier 
coefficients of x(t) at positive frequencies, while causing the 
coefficients at negative frequencies to vanish. 

A set of N unique amplitudes a n (t) and phases <fi n (t) is 
then implicitly defined by 

ai (t)e*« 



-(*) = 



"24xi(i)~ 




24x 2 (t) 




_2Ax N (t)_ 





a 2 (t)< 



i<fa(t) 



a N (t)e 



i4> N (t) 



(6) 



with the amplitudes being non-negative, a n (t) > 0. The 
nth amplitude a n (t) and phase 4> n (t) constructed in this 
manner are called the canonical amplitude and phase asso- 
ciated with the nth signal channel x n (t). Taking the real part, 
x(t) = 5i{x + (t)}, each signal channel is now described as a 
modulated oscillation with time-varying amplitude a n (t) and 
phase 4> n (t). The derivative of the nth phase, w n (t) = 4>' n (t), 
is called the nth instantaneous frequency fl], 0, ll22l . which 
gives the local frequency of oscillation of the nth signal. 

The analytic signal method provides the foundation for 
describing each channel of x(t) as an oscillation with time- 
varying properties. While the assignment of an amplitude and 
a phase to a given real-valued signal cannot be unique, the 
compelling properties of the amplitude and phase derived from 
the analytic signal are now well known |l|-[4], (6); see J5) 
for a useful review. Since a wide variety of physical processes 
can be aptly described as a set of modulated oscillations, the 
representation of the multivariate signal x(i) as in © is a 
strongly motivated and powerful model. 

III. Modulated Multivariate Oscillations 

In this section the notion of a modulated multivariate 
oscillation is formalized. The key is a local expansion of the 
signal about a demodulated version of itself. This expansion 
quantifies the signal's departure, at each moment, from the best 
possible fit to a set of sinusoidal oscillations all sharing a single 
frequency, i.e. from a pure oscillation. First-order and second- 
order deviations are introduced which quantify instantaneous 
linear and quadratic modulation, and which play a central role 
in an aggregate description of the signal's variability. 



IEEE TRANSACTIONS ON SIGNAL PROCESSING. SUBMITTED JANUARY 19, 2013 



Float Trajectories 



Signal Ellipses 



Residuals 



Bias Ellipses 





(b) 




s 

0' 



(d) 



27.5 25 22.5 20 
West Longitude 



30 27.5 25 22.5 20 
West Longitude 



30 27.5 25 22.5 20 30 27.5 25 22.5 20 
West Longitude West Longitude 



Fig. 1. Application of the extraction algorithm for modulated multivariate oscillations to freely-drifting oceanographic instruments from the northeast 
subtropical Atlantic 1161 . 1171 . The observed data x (t) in (a) is decomposed into a set of estimated modulated oscillations x(i), shown in (b) as a set of 
time-varying ellipses, plus a residual x r (i) in (c). In (b), ellipses are shown at twice actual size for presentational clarity. The time interval between the ellipse 
snapshots varies in time and is equal to the estimated instantaneous period, as described later in the text, with alternating grey and black ellipses. Panel (d) is 
the same as (c), but the ellipses represent the instantaneous estimated bias of the signal estimate. In (a) and (c), twenty-two different records are shown, with 
black lines used for those records for which a modulated oscillation is found, and grey lines for the remainder. The heavy gray curve in (a) and (c) outlines 
a particular record that will be used as an example later. 



A. A Local Signal Expansion 

The joint evolution of the multivariate signal x(t) in the 
vicinity of time t may be locally represented in terms of a se- 
ries of deviations from a set of constant-amplitude oscillations 
all evolving with some common instantaneous frequency uj(t). 
With r representing a time offset or "local time", the analytic 
signal may be written in the vicinity of a reference time t as 

x+(i + r) = e iw ^ T {x+(t) + rxi(t; w(t)) 

+ ir 2 x 2 (t; w(t))) + e 3 (t, t; w(t)))| (7) 

a representation we refer to as the local modulation expansion. 
The local modulation expansion describes the evolution of a 
multivariate signal as being due to the phase progression at a 
single time-varying frequency u(t), together with a series of 
deviations from this behavior. This model of joint structure is 
a key contribution, since it represents the multivariate signal 
as a single object, rather than as a set of unrelated oscillations. 

The p th vector- valued coefficient of the expansion, termed 
the p th-order deviation vector, is given by 



r)P r 



T = 



while the remainder term takes the form 

-* u W T x + (t + T) 



1 d 3 

ex (t,r;^))^-T 3 — 



(8) 



(9) 



for some (unknown) point u contained in the interval [0, r], 
as follows from the Lagrange form of the remainder in the 
Taylor series ll23l p 880]. To derive (|7), write 



x+(t + r) = e"'* 



*+(*• 



(10) 



and then Taylor-expand the term in square brackets with 
respect to the point r = 0. 

The local modulation expansion (0 states that the lowest- 
order joint behavior of the signal x+(i + t), considered as 
a function of local time r in the vicinity of a fixed reference 
time t, is for all signal channels to undergo a phase progression 
at a single frequency u>(t). The next-order behavior is a 
linear tendency in r, controlled by the first deviation vector 
xi(t;o;(t)), which is also subject to the phase progression at 
frequency u)(t). The still next-order behavior is a quadratic 
tendency in r, controlled by the second deviation vector 
x 2 (i;w(i)). Since xi(f;cj(i)) and x 2 (t;w(i)) are complex- 
valued in general, they impact both the amplitudes and phases 
of the oscillations in the various signal channels. In the vicinity 
of times t for which the signal x(i) is usefully described as 
a modulated oscillation at frequency u)(t), the remainder term 
e x (t,T;ui(t)) is expected to be negligible provided r is not 
too large compared with the oscillation period 2n/u>(t). 
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B. The Best Fit Frequency 

Our aim is to describe the common or joint oscillatory 
structure of the signal x(t). To this end, note that the quantity 



«*(t;w(t)) 



= \\^{t-Mt))\\ 

~ l|x + (t)|| 2 



_{t) - iu](t)x + (t)\ 
l|x+(i)|| 2 



(11) 



is the normalized instantaneous error involved in locally 
approximating the rate of change of the analytic version of 
x(i) as undergoing a uniform phase progression with some 
local frequency ui(t). For example, with x+(t) = a a e lulot , 
x' + (t) — iu x+(t) vanishes. One way to determine the best 
choice of frequency in the local modulation expansion (O 
is therefore to find that u)(t) which minimizes the error 
v x (t;uj(t)). Differentiating ( fTTT i with respect to ui(t) at each 
time t gives 



d 



2 d[u)(t) 



■«i(t;w(i)) 



l|x + W|| 2 



+ u(jt) (12) 



and we see, upon setting this quantity equal to zero, that an 
extremum in the fractional error v x (t;uj(t)) occurs for 



_ 9{xff(t)<(t)} 



(13) 



which is the power-weighted average of the N component 
frequencies 0J n (i). The second derivative of (fTTl i is positive at 
this value of u)(t), so this extremum is in fact a minimum. 

Thus Lu x (t) defined in ( fT3l minimizes the leading-order de- 
viation vector in the local modulation expansion (0, and is in a 
sense the "best fit" local frequency. The expression ([13} has in 
fact been encountered before, in [18|. Therein it was shown 
that the power- weighted time average of w x (t) satisfies an 
important global constraint — it recovers the first moment of the 
channel-averaged Fourier spectrum of x + (t) — and thus w x (i) 
generalizes the concept of "instantaneous frequency" (TJ, @, 
ll22l to the multivariate case. That this joint instantaneous 
frequency u) x (t) also has a compelling local interpretation 
as the solution to a minimization problem is another reason 
why it is a natural measure of the common time-varying fre- 
quency content of x(i). The interpretation of the instantaneous 
frequency as the solution to a local minimization problem 
holds for the standard univariate instantaneous frequency, since 
cu x (t) = cj)' x (t) for x+(t) = Oj^e^W is the TV = 1 special 
case of the joint instantaneous frequency. 

Henceforth we choose w(t) in the local modulation expan- 
sion (0 to take the value uj(t) — w x (t), that is, we write 



x+(t + r) = e 



iuj x (t)r , 



x + («)+rx 1 (t) + -T 2 x 2 (t) + e x (i,r) } (14) 



where the deviation vectors and residual are defined as 

x p (t)=x p (t;Lo x (t)) (15) 
e x (t,r) = e x (t,T;uj x (t)). (16) 



We refer to the x p (i) as the intrinsic deviation vectors, since 
the demodulation can be seen as a sort of coordinate trans- 
formation, with the natural or intrinsic choice of coordinate 
system being the one in which the phase progression follows 
the joint instantaneous frequency. 

The first two intrinsic deviation vectors are central in 
understanding the time-dependent joint structure of x(t) as 
a modulated oscillation. These are explicitly given by 

xi(t)=x' + (t)-iwx(t)x + (t) (17) 
x 2 (t) = x'|(t) - *2ufe(tX(t) - w 2 (t)x+(t) (18) 

the right-hand sides of which are oscillator equations that 
describe the first-order and second-order departure, respec- 
tively, of the evolution of x+(i) from a local oscillation at 
the frequency w x (t). The magnitudes these vectors, compared 
with the signal strength, are quantified by 



t(t)s llxtWH 



|x+W|| 



\\Mt)\\ 

|x + (t)|| 



(19) 



which will occur frequently in what follows. The first of these 
was also encountered in ifTSl , in which it was shown that 
v x (t) gives the time-varying contribution to the second central 
moment of the channel-averaged spectrum of x + (t) that is not 
explained by variations of the joint instantaneous frequency 
w x (t) about its time-mean value. Thus v x (t) is called the joint 
instantaneous bandwidth and is the natural multivariate gener- 
alization of the univariate instantaneous bandwidth introduced 
by Il24l - ll26ll . The local modulation expansion < TT~4T > shows that 
the joint instantaneous bandwidth v x (t) has a compelling local 
interpretation as the magnitude of the leading-order deviation 
of the multivariate signal x + (i) from oscillatory behavior. 



C. Physical Interpretation 

It is helpful at this point to say some words about the 
interpretation of the vectors that have been encountered. If 
x(t) is taken to represent a position, then x'(t) is a velocity 
and x"(i) is an acceleration, and x_|_(t), x! + (t), and x'!(t), 
are the analytic parts of the position, velocity, and acceleration 
vectors, respectively]! Then x x (t) could be termed the intrinsic 
analytic velocity, that is, that part of the analytic velocity 
which remains if the phase progression at the joint instan- 
taneous frequency is removed, and X2(i) could be termed 
the intrinsic analytic acceleration. It turns out that the third 
derivative of position, x"'(t), has an accepted name: it is called 
the jerk, see e.g. (27]. Thus the remainder e x (t, r) occurring 
in < TT~4T > is the supremum of the intrinsic analytic jerk. Con- 
straining e x (t, t) to be small therefore amounts to a kind of 
smoothness condition, namely, that the demodulated analytic 
signal does not exhibit too much jerkiness in its evolution. The 
fact that such smoothness is reasonable to expect for signals 
that may usefully be considered to be modulated oscillations is 
an argument in favor of our truncation of the local modulation 
expansion (fT4l at the quadratic term. 

2 Note that since differentiation and the analytic operator commute, the 
analytic part of a derivative is the same as the derivative of the analytic part. 
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D. The Deviation Vectors 

In this section we look at the deviation vectors in more 
detail. The squared norms of the deviation vectors take simple 
forms in the univariate case x+ (t) = a x (t)e i<fe W. Then flTJ 
for t! x (i) becomes \a' x (t)/a x (t)\, which is recognized as the 
modulus of the univariate instantaneous bandwidth [ 28 1 — [30 1 . 
Squaring ( fT8b for £ x (i) gives in the univariate case 



£(*) 



a x (t) 



+ K(t)} = 



a x (t) 



iu' x {t) 



(20) 



which involves a squared second derivative of both the ampli- 
tude a x (t) and the phase <fi x (t), since u x (t) = 4>' x (t). The 
complex-valued quantity a x {t) / a x {t) + iu' x (t), the squared 
magnitude of which occurs in d20t . has been previously 
identified as the coefficient of r 2 in the local modulation 
expansion of a univariate signal (9). A reasonable name for this 
univariate quantity is the instantaneous curvature. Then £ x (t) 
for N > 1 would be termed the joint instantaneous curvature, 
since it quantifies variability which has the same effect that 
amplitude curvature a x (t) and phase curvature 4> x (t) have in 
the univariate case. 

In what follows we will need some results concerning the 
deviation vectors. The first is that 



l|x+(t)|| 2 



-(t)W 



|x+(t)|| 



(21) 



and thus the quantity on the left-hand-side is real-valued. This 
states that the component of the first deviation vector along 
the direction specified by the signal vector is the fractional 
rate of change of the total signal amplitude. In the univariate 
case, this becomes the bandwidth a' x (t)/a x (t). To see ( f2~T~b . 
project the analytic signal onto its own first derivative to give 



x?(*x(f) i| x+ (i)ir 



iu x (t) 



(22) 



l|x + (t)|| 2 ||x + (i)| 
using the definition ( fT3l l of oo x (t) together with 

IM*)II' = |^xf (t)x+(t) = K{x* (*)<(*)} /||x + (t)|| 

and then (fJTJ follows from xi(t) = x' + (i) — iu) x (t)x + (t). 

The second result we will need is that £, x (t) > \uj' x (t)\, 
that is, that the instantaneous curvature is greater than the 
magnitude of the instantaneous chirp rate. To derive this, note 
that the derivative of the left-hand side of (l22l is 



d_ 

dt 



l|x+(i)|| 2 



x*(i)x"(i) , ||x' (t)|| 



|x+W|| 



|x+(t)|| 



- 2 



|x+(t)||'x^)x' + (i) 



(23) 



l|x+(i)|| l|x+(t)|| 2 
but we may also difference the right-hand side of (l22l to find 



d_ 

~dt 



r xf(t)xV(t) 
l|x+(t)|| 2 



d 
dt 



ix + wir 

l|x+(t)|| 



(24) 



The imaginary parts of these two expression combine to give 

^ /x?(tK(t)\ _ ... , o>- m i| X+ (t)ir 



|x+W|| 5 



) = ^ )+2WxW ™ (25) 



and using < TT~8T > to eliminate x" (t) together with (1221 1. we find 

•x£(t)x 2 (t)' 



|x+W|| 



(26) 



Now introducing the component of x 2 (i) projected onto the 
signal vector x + (t) as 

~ M xfft)x 2 (t) 

X2 '" W " |lx + (t)|| 2 X+W (2?) 
we may apply the Cauchy-Schwarz inequality, leading to 

|x 2 ,n(i)|| 2 > ||3{x 2 ,||(i)}|| 2 



(28) 



,|x + (t)|| 2 - ||x + (;)|| 2 

and hence £ x (t) > \uj' x (t)\, as stated. 

E. Pure Oscillations and Phase Signals 

To understand the distinction between the linear and 
quadratic terms in the local modulation expansion, we intro- 
duce two particularly simple types of multivariate oscillatory 
signals. A signal x(t) may be said to be a multivariate pure 
oscillation if its analytic part is given by 



x+(i) 



(29) 



for some fixed vector x Q and fixed frequency u> . Similarly 
x(i) may be termed a multivariate phase signal if 



x+(i) 



(30) 



for a fixed vector x Q and analytic phase function e 1 ^^. The 
multivariate phase signal d30b is the natural generalization of 
the univariate phase signal of e.g. 0. A univariate phase signal 
may be written as x + (t) — laole 1 ^*- 1 where the \a \ is the sig- 
nal amplitude and e 1 ^"^ is analytic. In the multivariate case, 
x is complex-valued in general as it incorporates information 
on phase shifts between channels. Thus the constant part x Q 
of the phase signal can only be interpreted as an amplitude for 
the univariate case N = 1, in which case it can be made real- 
valued and nonnegative by absorbing its phase into e l( ^ x (*' . 

Phase signals of the form (130b are frequency modulated, 
with each signal channel having identical time-varying in- 
stantaneous frequency u> n (t) = uj x (t) = <fi' x {t), but they are 
not amplitude modulated since all of the N amplitudes are 
constant. Phase signals are the more general class since all 
pure oscillations are phase signals but not vice-versa. Note 
that there are strong constraints on the class of phase functions 
</> x (£) such that e 1 ^^ be analytic; see e.g. the detailed 
discussion of univariate phase signals in Q. 

The intrinsic deviation vectors take very simple forms for 
these two types of signals. For a pure oscillation, x p (i) 
vanishes identically for all p > 0. For a phase signal, we have 



xi (t) 
*2 (t) 



«4(t)x+(t) 



(31) 
(32) 



so that the first deviation vector vanishes, but the second 
deviation vector is nonzero whenever the joint instantaneous 
frequency varies with time. The former expression (13~T1 follows 



IEEE TRANSACTIONS ON SIGNAL PROCESSING, SUBMITTED JANUARY 19, 2013 



7 



directly from (fTTT i. The latter (|32t may be readily found by 
rewriting JT8l for the second deviation vector as 



x 3 (i) = *«£(t)x+(t) + [x' x (i) - iw x (t)xi(t)] 



(33) 



where the first term depends on the joint chirp rate ui' x (t), and 
the other terms vanish when xi(i) vanishes. 

This illustrates a subtle distinction between the linear and 
quadratic terms in the local modulation expansion. Both pure 
oscillations and phase signals have vanishing linear deviations 
from local oscillatory behavior at all times, as measured by 
the norm of the first deviation vector xi(t). However, phase 
signals differ from pure oscillations at second order, since 
the former have non-vanishing quadratic deviations from local 
oscillatory behavior. Conversely, when Xj.(i) is negligible in 
the vicinity of time t, we may say that the signal locally 
evolves as if it were a phase signal having a frequency ui x (t). 
When X2(f) is also negligible, we may say that the signal 
behaves as a pure oscillation up to second order. When the 
leading-order term e lW:x; ^ T x+(t) dominates in (TBI for r not 
too large, we may say that the signal evolves in the vicinity 
of time t as would be expected for a pure oscillation having 
frequency w x (£). 

F. Definition of a Modulated Multivariate Oscillation 

We are now in a position to formalize what is meant by a 
modulated oscillation in an arbitrary number of dimensions. 
This is accomplished by proposing a single measure of the 
degree of departure of a multivariate signal from a pure 
oscillation. An TV-channel real-valued zero-mean signal x(i) 
is assumed to have an analytic version x+(t) that is defined 
and thrice differentiable over some time interval T. The 
analytic signal x+ (t) is then expanded via the local modulation 
expansion (fT~4T > using the joint instantaneous frequency ui x (t). 

Definition 1: The Modulated Multivariate Oscillation 
Let the local stability level 5t be the smallest positive constant 
satisfying for all t € T the constraints 









W x (t) 







together with 



sup ■ 



1 



M*,r)|| 



:|x+WII 



< Si 



< S T . 



(34) 



(35) 



ter |w x (i)r 

Strongly modulated signals correspond to large values of St, 
while St vanishes for a pure oscillation. The signal x(t) is said 
to be a modulated multivariate oscillation over time interval 
T if the local stability level is less than unity, St < 1. 

In this definition, modulated multivariate oscillations occupy 
a continuum — classified according the local stability level 
St — with pure oscillations as the limiting or ideal case of 
vanishing modulation strength. Signals for which St exceeds 
unity present temporal variability that locally exceeds the rate 
of change of phase at least somewhere on the time interval T. 
Such extremely strong modulation would be evidence that 
the signal is not well modeled in terms of an oscillation 
at a common time-varying frequency. An important point is 
that the class of modulated multivariate oscillations is far 



larger than that of the so-called "asymptotic" signals, see e.g. 
the discussion in Q, which roughly correspond to univariate 
signals having negligible modulation strength, St *C 1. 

In the next section, this ability to quantify the degree of 
variability becomes essential for determining the time-varying 
bias involved in the recovery a modulated oscillation from a 
noisy observation. 

IV. Multivariate wavelet ridge analysis 

In this section, a local optimization method — multivariate 
wavelet ridge analysis — is created that is able to extract esti- 
mates of a modulated multivariate oscillation from potentially 
noisy observations. Time-varying forms for the leading-order 
bias effects are also derived. This extends the work of [7 1 and 
[8 1 on the univariate wavelet ridge method, and that of [9| on 
its bias properties, to the multivariate case. 

A. Wavelet Basics 

To isolate a signal of interest from surrounding variability, 
a time/frequency localized filter is necessary. A wavelet ip(t) 
is a square-integrable complex-valued function satisfying the 
admissibility condition OD 

2 

■ duj <C oo (36) 



and the wavelet is said to be analytic if its Fourier transform 
^(uj) = J ip(t) e~ tult dt vanishes for all negative frequen- 
cies. The wavelet transform of a real-valued square-integrable 
vector-valued signal x(t) with respect to the wavelet tp(t) is 



"X,V> 



(t,s) 



1 ,» 

-ip 

s 



x(t) dr 



1 r°° 

— / y*{sLj)X(uj)e iut duj 
2tt Jo 



(37) 



(38) 



where the latter form follows by the convolution theorem. With 
the 1/s normalization in (f3Tb rather than the more common 
1/t/s, the wavelet transform can be seen as a set of bandpass 
operations indexed by the scale s. Considered as a function 
of time at each scale s, the wavelet transform is seen as a 
stack of analytic signals. The frequency-domain wavelet 
obtains a maximum modulus at a frequency uj, called the peak 
frequency. Without loss of generality, we set *§>(lo^) = 2. With 
these choices, the wavelet transform w x ^(t,s) of a sinusoid 
x(t) = <zi cos(wii 4- 4>\ ) obtains a maximum modulus at the 
scale s — and the value of this maximum recovers the 

amplitude of the sinusoid, \w Xt ^(t,uj^,/uJi)\ = \ai\. 

B. Joint Ridges of a Multivariate Signal 

The detection of a modulated oscillation within the analytic 
wavelet transform w X) ^(i, s) is accomplished as follows. A 
ridge point of w X) ^,(f, s) is a time/scale pair (t, s) satisfying 

— \\w^(t, s)\\ =0, — ||w x ,^(t, s)|| < 0. (39) 

Thus ridge points are locations where the norm of the wavelet 
transform vector achieves a local maximum with respect to 
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scale. This definition of a multivariate ridge poinj| is the 
natural generalization of the definition for the univariate case 
|9l , as proposed in the earlier preliminary work (TTOl . 

Adjacent ridge points are then connected to each other 
to yield a single-valued, continuous function of time called 
a ridge curve s~(i) that extends over some time interval T. 
In practice, two numerical thresholds must be introduced: a 
bound on the magnitude of to avoid jumps across scale, 

and a minimum ridge duration, to avoid spurious ridges that 
are short compared with the wavelet length. Having identified 
a ridge curve s"(t), the ridge-based estimate of the analytic 
signal is then given by 

x+(t) = w Xl ^(t,«(t)), teT (40) 

which is simply the set of values taken by the wavelet 
transform along the ridge curve. In order for the wavelet ridge 
estimate x+ (t) to be a good estimate, it is necessary to choose 
the wavelet properties to match the signal properties. This is 
addressed in the next section. 

An example of the multivariate wavelet ridge algorithm 
is shown in Fig. [2] One of the bivariate time series from 
Fig. [T] — specifically, that trajectory which is marked by the 
heavy gray curve in Fig. QJ and Fig. QJ — is presented to- 
gether with its wavelet transform using a choice of wavelet 
and parameter settings to be discussed later in Section fVl 
The wavelet transform modulus ||w Xl i/,(t, s)|| shows a clear 
maximum value as a function of scale, expressed here as 
the period 2tts/uj^,. The scale at which this maximum value 
occurs changes considerably throughout the record, decreasing 
by an order of magnitude from the beginning to the end 
as well as presenting some low-frequency variability. The 
multivariate ridge curve s(t) is seen to follow the variability of 
the maximum of ||w Xj ^,(t, s) || as a function of time. Evaluating 
the wavelet transform along this time-dependent curve as in 
( f4Qb defines the estimated modulated oscillation x+(t), which 
is plotted in Fig. Q] as a set of time-varying ellipses together 
with the estimated residual x r (t) = x Q (t) — 5R{x + (t)}. 

In the following we will use a measure of the distance, at 
time t, of a scale point s from the instantaneous frequency 
curve ui x (t), called the scale deviation 

Aw x ^(t,«) = — ^-1. (41) 

On the time-varying scale curve s(t) = u^,/u x (t) corre- 
sponding to the instantaneous frequency curve w x (i), the 
scale deviation vanishes. One may envision that the constraint 
|Aw x . ,/,(£, s)\ < \c\ for some small \c\ > paints out a swath 
surrounding the instantaneous frequency curve, with the width 
of this swath increasing as \c\ increases. If |Ao> Xi w,(t, s)| < 6%,, 
the scale point s is said to lie in the neighborhood of the 
instantaneous frequency curve at time t. In the next sections we 
find the conditions under which the ridge equations ( |39l > have 
a solution within the instantaneous frequency neighborhood. 

3 This type of ridge point is called an "amplitude ridge point" by |9|. 
Another possibility is a "phase ridge point", utilizing a stationary phase 
condition as proposed by |8 j. However, since \9\ find negligible difference 
between the two types of ridge points in a perturbation analysis, and give 
practical reasons to prefer the amplitude ridge points, only these will be 
considered here. 



Multivariate Ridge Method Example 




100 150 200 250 
Day of Year 1986 

Fig. 2. A bivariate position signal, differentiated in time for presentational 
clarity, is plotted in (a). The solid curve represent eastward velocity and 
the dashed curve represents northward velocity. The modulus of the wavelet 
transform w x ^(t, s) shown in (b), has units of kilometers and is plotted with 
a logarithmic y-axis. The contours range from to 65 km with a spacing 
of 5 km. The heavy curve is a single unbroken ridge resulting from the 
application of the multivariate ridge algorithm. This time series is from the 
position signal marked by the heavy gray curve in Fig. [TJ and Fig. [TJ. 



C. Constraints on the Wavelet 

The most important wavelet parameter after its peak fre- 
quency u^p is the dimensionless duration P^, defined by 



P,l, = 4 / -U) 



(42) 



The quantity under the radical is positive for a wavelet with 
a real-valued Fourier transform ^(u>), since the wavelet then 
obtains a maximum value at uty, making ^"(uty) negative. 
It may be shown that P^/tt corresponds to the number of 
oscillations at period that fit within the central window 

of the time-domain wavelet ll32l . Also is seen as a 

dimensionless measure of the wavelet bandwidth, since a 
Taylor expansion of &(u)) about the peak frequency gives 

2" 



1-i 

2 



'■UtP- (43) 



If Pjf, — 1, the half-power points in this quadratic approxi- 
mation occur at zero and 2uty, and so the frequency support 
is extremely broad. One thus expects P^ > 1 for wavelet 
functions that are usefully localized in the frequency domain. 

More generally, the basic features of the wavelet may be 
characterized by its dimensionless derivatives B32II 



(44) 



which are then evaluated at the peak frequency ojj,. Note that 
^(uty) vanishes by definition and that P^ — — ^(w^). The 



wavelet suitability criteria Q9J 



1 

pi 



< 



-p/2 

-(P-1J/2 



fez 
p+i 

2 



€ Z 



(45) 



are a set of conditions that limit the size of the wavelet's 
dimensionless derivatives at the peak frequency u]$, compared 
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with the inverse of the signal stability level St- As the signal 
becomes more rapidly varying, St increases, and these bounds 
on the size of the dimensionless derivatives becomes more 
stringent. Note that the second condition in (PES! places a 
stronger constraint on odd derivatives compared to the even 
derivatives, a reasonable constraint that also proves convenient 
for the subsequent analysis; see [9| for details on this choice. 

It will be seen shortly that the wavelet suitability conditions 
are the key to ensuring that the bias of the ridge-based signal 
estimate remains small. The lowest-order suitability condition, 
at p = 2, implies the wavelet duration is P^, < ^2/ St- 
This means that the more strongly the multivariate signal x(t) 
varies over an oscillation period — reflected by an increasing 
value of St — the fewer oscillations the wavelet can contain 
in time. But as the wavelet becomes narrower in time it must 
become broader in frequency, and combining this with the 
earlier discussion, one expects 1 < P^, < y/2/Sr for wavelets 
that are usefully localized in both domains. If St were to 
become large compared to unity, no such wavelet could be 
found. Our definition of a modulated multivariate oscillation 
implies St < 1, ensuring that the suitability conditions may 
be satisfied for a wavelet with P^ in the range 1 < P^ < 

D. The Wavelet Transform of a Multivariate Oscillation 

The wavelet transform w Xj ,/,(i, s) can be expanded into 
a hierarchy of terms that reveal the interaction between the 
signal variability and the wavelet shape, allowing the lowest- 
order bias term to be identified. A change of variables ap- 
plied to the wavelet transform (137) . and substituting x(t) = 
[x + (t) + x* + (i)] /2, leads to the form 



1 f°° 

w x ,v>(*, s) = - / ^*(T)x + (t + sr)dT 

* J-oo 



(46) 



where the contribution of x*^ (t + r) to the integrand vanishes 
on account of the analyticity of the wavelet, as is clear form 
the Fourier-domain form ((38). Inserting the local modulation 
expansion (fl4] > of x + (t), we obtain 

^ J-oo 

c+ (t) + st x i (t) + i (sr) 2 x 2 (t) I dr 

+ Aw x ,^(i,s). (47) 

Note that to guarantee the square-integrability of the r 2 term in 
the modulation expansion, the long-time decay of the wavelet 
must be |^(t)|/|^>(0)| ~ \t\~ r '*• for some number > 3. 

There are some subtleties surrounding the residual term 
Aw x ^(t, s). This is implicitly defined as the difference 
between the left-hand side of d47l i and the integral on the 
right-hand side. It is not the same as the wavelet transform 
of the Taylor-series remainder e x (i, r), because the form (0 
for e x (t,r) is only valid over the interval T where we have 
assumed the signal is differentiable. Bounding the residual 
term Aw x ^ (t, s) has been examined in the univariate case by 
(9), and since there is no major difference in the multivariate 
case, we refer the reader there for a detailed discussion. In 



general we may expect this term to be very small when the 
time decay of the wavelet is stronger than t~ 3 for signals that 
meet our definition of modulated multivariate oscillations. 

Assuming the suitability criteria (l45l l are satisfied, the 
wavelet transform in the instantaneous frequency neighbor- 
hood, Aaj x .^(i. s) — 0(S T ), takes the simple form 

0(5 T ) 



X 2 (*) 



0(& 2 T ) 



w x ,0(M) = x+(t) - a 



(t) 



Aw x ,^(t,«) (48) 



as will be proved shortly. This powerful result states that 
at scale points sufficiently close to instantaneous frequency 
curve, the wavelet transform approximately recovers the ana- 
lytic signal x + (i). The time-varying forms of the two lowest- 
order deviations from the analytic signal, up to second order 
in St, are explicitly resolved. 

The derivation of (l48l is as follows. Note that the pth 
frequency-domain derivative of the wavelet is given by 



*W(t, 



(-ir) P V(r)e" 



r dr. 



(49) 



Evaluating the conjugate of this quantity along the time- and 
scale-varying frequency cj = scu x (t), we may define a joint 
function of the wavelet and the signal as 



(irf ^*(r)e 



isuj x (t)r 



dr (50) 



which is a function of the time-scale plane. Inserting 
the wavelet transform expression (l47t leads to 



into 



w x ,^(t, s) = $ (t, s)x + (t) - i$i(t, s) 



gift) 



-i$ 2 (t, s )^+Aw x> ^(t, S ) (51) 
2 w x (i) 

in which the first two deviation vectors appear explicitly. 
This can be simplified by finding approximate expressions 
for the Q p (t, s) that are valid in the instantaneous frequency 
neighborhood. 

In terms of its dimensionless derivatives ^ p (oj), the wavelet 
has a Taylor expansion about the peak frequency ui^ of 



°° 1 ~ ( s 



fe=0 



SU) 



(52) 



Differentiating both sides, and recalling &(uj^) = 2, we find 



(53) 



IEEE TRANSACTIONS ON SIGNAL PROCESSING, SUBMITTED JANUARY 19, 2013 



10 



for the Taylor series expansion of the pih derivative of the 
wavelet, after employing a change in the index of summationQ 
Thus d50l l becomes, after making use of (|4TT >. 

$ P (M) = [l + A^(t,s)fx 

GO - 

^^ +p K)[A^(M)f (54) 

k=0 

in terms of the scale deviation. This may be approximated 
in the instantaneous frequency neighborhood by making 
use of the wavelet suitability conditions, and recalling that 
Aui Xt ^(t,s) is 0{8 T ) in the instantaneous frequency neigh- 
borhood by definition. The first five $ p (t, s) are found to be 



$ (M) 


= i + o(4) = 0(1) 


(55) 


*i(t,a) 


= Acj x ^(t, «)*S(uty) + 0(4) = 0{5 T ) 


(56) 




= %(^) + 0(8 t ) = 0{5t 1 ) 


(57) 


$s(*,«) 


= %{^) + 0{l) = 0{5- 1 ) 


(58) 


$ 4 (M) 


= ^K^p) + 0(S T ) = (<5 T 2 ) 


(59) 



in the instantaneous frequency neighborhood. Using (I55l-(l57l 
and gathering terms by order in (IBTi . the result (l48l follows. 

E. The Bias of the Wavelet Ridge Method 

One may now show that the ridge equations have a solution 
within the instantaneous frequency neighborhood. That is, 
there exists a scale curve 



u x (t) 



[i + o(4)] 



(60) 



which satisfies the ridge equations (139b . The proof of this 
statement, which extends a similar result for univariate case 
by J9] to the multivariate case, is given in Appendix IDl Then 
within the instantaneous frequency neighborhood, (|4F1 applies, 
and we obtain 



(f) 



,V, (i,s(t)) = x+(i) 



1 D 2 X 2 (Q 



rP, 



+ 0(4) + Aw x> (t, w^/w* (*)) (61) 



as an explicit form for the estimated signal x^(i), resolving 
the lowest-order time-dependent bias term. 

There are several important implications of this result. The 
leading error term is due to the second deviation vector, 
and not the first deviation vector. This is attractive since it 
means that the lowest-order deviation of the signal from a 
modulated oscillation — a linear tendency in local time r — 
does not impact the analysis errors at lowest order. This 
leading-order error is seen to be associated with the value 
of the wavelet transform along the instantaneous frequency 
curve, with the deviation of the ridge from the instantaneous 
frequency curve only contributing at higher orders in St- A 
measure of the total error of the signal estimate is the norm 

4 Note that {53) corrects the similar expression (1 14) of [9 . The latter is only 
approximately correct, up to order 5~ . Subsequent perturbation expansions in 
|9| are not affected because the erroneous contributions occur at unresolved 
orders in Sj<- 



of the difference between the original signal and the estimate 
normalized by the signal amplitude, found to be 

||x V; (t)-X+(i)|| _ lrflfixWI 



2 



(62) 



|x + (i)|| - 2-*a£(t) 

which is controlled by the joint instantaneous curvature. Since 
Pjf, is the wavelet duration, this states that the error is 
proportional to the degree of signal curvature over the time 
support of the wavelet. To make the leading-error bias term 
negligible, one must be able to choose the wavelet duration 
such that this term is sufficiently small. 

The joint instantaneous frequency cu x (t) may also be esti- 
mated. One possibility is to form an estimate by substituting 
the signal estimate x^, (t) for the true signal x + (t) in the defi- 
nition of the instantaneous frequency ( fT3l . However, following 
|9| , we instead form the joint transform frequency 



{<(M)|w x ^(t,s)} 



(63) 



||w x ,v(M)|| 

and evaluate this quantity along the ridge to obtain the ridge- 
based instantaneous frequency estimate 

Q x ^{t) = n x ^(t, s(t)). (64) 

In Appendix IDl we find 



1-rP, 



1 u2 3 {xfft)x 2 (t)-xg(t)x 3 ft)} " 

||x+(t)M(t) 
= io x (t) [1 + 0(4)] (65) 



as an expression for the time-dependent form of this instan- 
taneous frequency estimate. The wavelet suitability conditions 
ensure that this estimate is accurate to second order in the 
local stability level St- 

The bias itself may be similarly estimated. We form the a 
version of second deviation vector associated with the wavelet 
transform 



d 2 

W 2 ;x,v(M) = ^2 W x,«HM) 



d 

i2Q^(t,s)—w x ^(t,s) 

-^(*,»)wM(i,«) (66) 

which is created by substituting w x ^(t, s) for x+(£) in ( fl~8b . 
replacing total time derivatives with partial time derivatives. 
Then we have the estimates 



X2;V>(*) = W 2; x,i/>(*>s(*))j §t,v(*) 



X 2;V(*) 



(67) 



llx^WH 

for the second deviation vector and its modulus, the joint 
instantaneous curvature. This permits the bias of the estimated 
signal, and the normalized bias magnitude in j62l , to be 
estimated. 

V. Application 

This section illustrates the multivariate wavelet ridge 
method with an application to real-world bivariate data. The 
data is from a set of instruments tracking the ocean currents, 
and is representative of a large amount of similar oceano- 
graphic data, see e.g. ifTTl . 
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A. Data 

The data, shown earlier in Fig. QJ, consists of position 
records from twenty-two freely-drifting acoustically-tracked 
subsurface floats. These were deployed off the west coast of 
Africa in the eastern North Atlantic in order to study the local 
currents in an early experiments of this type [1161 . IfPTl . The 
instruments are designed to remain neutrally buoyant near a 
particular depth, 1000 meters in this case, and are tracked 
acoustically by triangulating sound travel times between the 
instruments and nearby fixed points. In the experiment shown 
here, the sample rate was one day, and only float records with 
a length exceeding 200 days are presented. This dataset and 
many similar ones are available from the World Ocean Circu- 
lation Experiment Subsurface Float Data Assembly Centerjf] 

B. Choice of Wavelet Family 

In implementing the wavelet ridge analysis, the choice of 
family of analytic wavelets emerges as being important to 
obtaining desirable properties of the transform, an issue that 
has been investigated in detail by |32| and [9|. A particularly 
attractive choice is the generalized Morse wavelet family [32|- 
||34l , given by the frequency-domain from 

*/3, 7 M = J7(w)a /3 , 7 w /3 e-" T (68) 

where U(ui) is again the unit step function, ap n is a normal- 
izing constant, and B and 7 are two adjustable parameters. 
The peak frequency occurs at Lup n = (/3/7) 1 / 7 , and we 
choose ctp y = 2(e7//3)^/ 7 in order to meet the convention 
^ p n {uj p :1 ) = 2. For the generalized Morse wavelets, the 
duration takes the simple form Pp a = \/^7- in [ 32 1, the 
7 = 3 family is recommended as a superior alternative to 
the only approximately analytic Morlet wavelet. With this 
choice, the wavelet duration Pg, 7 is matched to the signal 
variability by adjusting B. It is also shown in l32| that time 
decay of the generalized Morse wavelets is controlled by B, 
with 1^,7(^)1/1^/3,7(0)1 ~ Thus to ensure square 

integrability in ( f4Tb . the constraint > 3 translates to (3 > 2. 

C. Ridge Application 

The multivariate wavelet ridge analysis method using the 
generalized Morse wavelets is applied to the data set shown in 
Fig- QJ, using a freely distributed software package described 
in Appendix lAl For all but two of the time series, the 7 = 3, 
fi = 3 generalized Morse wavelets are used, so in this 
case we have Pp n — y/]3j — 3. The wavelet transform 
vector w x . s) is computed with 82 logarithmically spaced 
frequency levels, with a lowest frequency of 0.01 cpd (cycles 
per day) and a maximum frequency of 0.28 cpd. For two of the 
time series, the frequency content was at considerably higher 
frequencies, and so we use different settings in order to more 
closely approach the Nyquist frequency. For these two time 
series, we used the 7 = 3, 8 = 8 generalized Morse wavelets, 
so Pp n — 2V6, and computed the wavelet transform at 140 
logarithmically spaced levels with a lowest frequency of 0.01 
cpd and a maximum frequency of 0.34 cpd. 

; http://wfdac.whoi.edu 



The multivariate ridge method described in Section II V-B I 
is then applied, rejecting ridges with that execute a smaller 
number of complete cycles than 2Ps ~ = 6. At a very small 
number of points, two valid ridges are obtained which overlap, 
and these are combined into a single estimated signal x^(i) 
through a power-weighted average. Thus there is either one or 
zero estimated modulated oscillations x^(i) present at each 
time. These modulated bivariate oscillations can be converted 
into the parameters of a time-varying ellipse following [18|, 
and these ellipses are shown in Fig. QJ. The time interval 
between successive ellipses is proportional to the estimated 
period 2ir/u) x (t), and the ellipses are shown at twice actual 
size. A set of estimated residuals formed by subtraction, 
x r (t) = x Q (i) — 5i{x,/,(i)}, shown in Fig. Finally, the 
estimate bias is shown in Fig. [TJl by converting the estimated 
deviation |Pjx2 ; ^(i) into time-varying ellipse parameters. 
That the estimated bias is generally small compared to the es- 
timated signals is consistent with visual inspection of Fig.Q]:, 
in which the residual curves appear to be largely devoid of 
oscillatory motions. 

VI. Conclusion 

This paper has addressed the analysis of modulated oscil- 
lations in multivariate time series. The key contribution is a 
local expansion of modulated oscillatory variability in terms 
of deviations from a pure oscillation at a common but time- 
varying frequency. This model captures the essence of time- 
dependent wavelike motion spanning multiple signal channels. 
A condition for a signal to be considered a modulated multi- 
variate oscillation is given, which amounts to demanding that 
the magnitude of the local deviation of the signal from a pure 
oscillation is no larger than the magnitude of the signal itself. 

A generalization of wavelet ridge analysis for multivariate 
timeseries is presented which enables an estimate of the 
modulated oscillation to be formed from a wavelet transform 
of the signal. By appealing to the signal model, constraints 
may be placed on the choice of analyzing wavelet such that 
the estimate of a modulated oscillation is guaranteed to have 
small bias. By considering signals which are both multivariate 
as well as non-negligibly modulated, and by presenting forms 
for an important source of time-varying error, this work 
substantially extends earlier tools for analysis of nonstationary 
or modulated oscillations. 

Appendix A 
A Freely Distributed Software Package 

All software associated with this paper is distributed as a 
part of a freely available Matlab toolbox called Jlab, written 
by the first author and available at http://www.jmlilly.net 
The Jsignal module of Jlab includes numerous routines for 
multivariate wavelet ridge analysis suitable for large data sets. 
The wavelet transform using generalized Morse wavelets is 
implemented with wavetrans, which calls morsewave 
to compute the wavelets. The standard univariate and joint 
wavelet ridges are found by ridgewalk using a numer- 
ically efficient algorithm that includes quadratic interpola- 
tion between discrete scale levels. Position records given 
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in latitude and longitude are converted into displacement 
velocities with latlon2uv, while and latlon2xy and 
xy21atlon convert between latitude and longitude and a 
local Cartesian coordinate system. Ellipse parameters are 
found from a pair analytic of signals with ellparams, and 
the ellipses can be plotted using ellipseplot. Finally, 
makef igs_multivariate generates all figures in this 
paper. 

Appendix B 

Expansion of the Rate of Change of the Signal 

In this appendix a version of the local modulation expansion 
([T4l > for the first time derivative of the signal is derived. The 
partial derivative with respect to the global time t of x+(f +t), 
assumed thrice differentiable, may be expanded as 



d_ 

dt 



*+(*- 



+ ir 2 [x 3 (t) + iuj x (t) X2 {t)} + i-c4(t)T 3 x 2 (i) 



+£x ;t («,r)} (69) 

where e x;i (t,r) is a remainder term. To derive this, write the 
partial i-derivative of x + (t + r) as 



d s O r 



d 



3 i(jj x (t)r 



,-«W x (t) 



T x+(t + r)]} 



dt 



e -«-x(*)r x+ (f + T )l + i TU ' x (t) x+(t + r). (70) 



Substituting from (fl4] i. the term in square brackets on the 
second line can be Taylor-expanded in r as 



d_ 

dt 



«(*)•> 



x+(£ + r) 



= x' + (<) + ™i(t) + -r 2 Z 2 (t) + e x ,(t,r) (71) 



where the residual takes the form [23. p 880] 



e x '(*> T ) 



d 3 f 9 r 



1 



6 <9t 3 at 



-2LJ x (i)r 



x+(t- 



(72) 



for some (unknown) point v contained in the interval [0, r]. 
Note that v is not in general the same as the point u 
appearing in the expression (0 for the remainder e x (t, r) in 
the comparable expansion of x_|_(t). The identities 

x' x (i) + «4(i)x+(i) = x 2 (t) + ^w x (^)X 1 (^) (73) 
X 2 (i) + 2iwi(t)xi(t) - x 3 (t) + iwx(t)x a (t) (74) 

may readily be verified from the definitions ((8) of the deviation 
vectors. Substituting these into (|7T1 i. and combining the result 
into d70l > together with the local expansion in r of x + (t + r) 
given by (fl4l . we obtain d69l with 



(75) 



e x; t(t,T) = e x /(i,r) + iruj' x (t)e x (t, r) 



as the form of the remainder term. In the above, we have taken 
care to avoid differentiating a remainder term such as e x (t, r). 



Appendix C 
Time Derivative of the Wavelet Transform 

Here an expression for the time derivatives of the wavelet 
transform is found that is valid near the instantaneous fre- 
quency curve, using the expansion d69l for the rate of change 
of the signal derived in the previous appendix. The normalized 
time derivative of the wavelet transform is found to be 



(t,s) 



1 1 



c(*) 



c (t)2 



r(r)^ + (t 



x' (t) 

$o(i,s)^777-^i(M) 



U) x (t) 

x 3 (t) ,.x 2 (t) 



f st) dr 

Xi(i)' 



2 <J>3 ^' Sj a; x (t) Wx (t) 
+ Aw xrf;( (M) (76) 



by inserting d69l into the time derivative of (|46] i, exchanging 
the orders of differentiation and integration, and making use of 
the definition d50t of the $ p (t, s) functions. The residual term 
Aw^^tfi, s) here is again implicitly defined as the difference 
between the left-hand side and the other terms on the right- 
hand side; see Appendix D of [9| for details on bounding this 
term. Gathering orders in (l76l . we find 




xiW 
+ 0(5 3 T , 



Aw xi;( ((,s) (77) 



making use (T55l>— d58l> together with fact that c<J x (i)/c<J x (t) is 
0(6%,)- This implies that the first derivative of the signal 
x^ (i) is accurately estimated by the value of the partial time 
derivative of w x ^(i, s) evaluated along the ridge curve. 

For what follows, assume that the wavelets are real-valued 
in the frequency domain. Note that in the instantaneous 
frequency neighborhood the imaginary part of the projection 
of the wavelet transform onto its own time derivative is 



1 



w x (t) 



l!x+«ll 2 

0(S T ) 



~ 1 S{xf(t)x 2 (f)} 1~ g(t) 

o(4) 



9f{xf(t)x a (t)-x?(t)SE3(t)} 



+ o(4) + o 



l|x + (t)|| 2 



l|x+ 
Aw 



|x + («)|| 



(78) 



as follows by combining the two expansions (|48T > and ( fTTT i, and 
using the substitution Xj_(t) = x+(t) + ioj x (i)x + (i); note that 
the order <5t term in ( l78l arises twice, and thus its coefficient 



IEEE TRANSACTIONS ON SIGNAL PROCESSING, SUBMITTED JANUARY 19, 2013 



13 



is unity rather than 1/2. At the same time we may find, again 
for real-valued wavelets, 

O(St) 



9?{x£(i)x 2 (i)} 



-(*)ll 



o(4) 



1 



£(*) 



-(*)« 



Aw x ^(t,s) 



for the expansion of the modulus-squared wavelet transform 
in the instantaneous frequency neighborhood. In deriving both 
of these expressions we have made use of the fact that 
x^(£)xi(i) is purely real, see (|2"TV 

The transform instantaneous frequency in the <5^ neighbor- 
hood of the signal's instantaneous frequency is then given by 



w x (i) w x (t) 



|w x ^(t,s)|| 2 
o(4) 



1 + 5* a ^^ 



{xf(t)x 2 (t)-x*(i)x 3 (i)} 



0(4) + o 



Aw x ^(i,s) 



+(*)« 
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(80) 



| x+ (t)ip ; v llx+WH 2 

which is found by combining (F78b and ( |79l l, using 1/(1 + a;) = 
1 — x + x 2 + 0(x 3 ), and carefully keeping track of the cross 
terms from the product of the two expansions. A number of 
cancelations occur: the leading-order term in the numerator 
cancels the leading-order term in the denominator, and the 
square of the first term in ((79) (arising from the expansion of 
the denominator) cancels an identical term arising from the 
product of the numerator with the expansion of the denom- 
inator. The result (180) implies that the transform frequency 
evaluated along the ridge, u)^(t) = H.^(t, s(i)), is an accurate 
estimate of the joint instantaneous frequency w x (t). 

Appendix D 

Scale Derivative of the Wavelet Transform 

The scale derivative of the wavelet transform can be found 
in a similar fashion to the time derivative in the preceding 
appendix. The scale derivative of the shifted analytic signal 
x + (t + st) is related to its time derivative via 



s— x+ (t 

os 



9 , 
ST dl X+{t 



Then inserting 
transform expression 



st). (81) 
into the scale derivative of the wavelet 



and again using (IBTTt leads to 

d 1 f°° d 

s^w XlV ,(t,s) = - J st%1>*{t) —x + (t + sr)dT 



x' (t) 
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W 3(t) W»(f) 
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<(t) x 2 (t) 
a£(t) w a (t) 
+ Aw xAs (M) (82) 



with the residual Aw Xi w,. s (t,s) defined implicitly as before. 
Again, we refer the reader to Appendix D of |]9] for details 
on bounding this term. Noting (|55ll-(l59l and using x' + (t) = 
xi(t) + iuj x (t)x + (t), we can gather orders in (l8Zt to yield 



O(l) 



s T w x ^t,s = -«* 2 (w. — 

os uj x (t) 
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x 2 W 
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*2K) + 2*sK) 

+ 0(4) + Aw^(t,«) (83) 

in which terms up to first order in St are resolved. This implies 
that a suitably normalized version of the scale derivative of the 
wavelet transform evaluated along the ridge recovers the first 
intrinsic deviation vector Xi(i). 

The ridge condition can now be evaluated using expressions 
for the wavelet transform and its scale derivative. The ridge 



condition ^||w Xi ^,(i, s)\\ — 0, from ( [39l , is equivalent to 



R \ llx+WII 2 
and inserting (l48l i and d83l , we find 
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(85) 



Assuming that the wavelets are real-valued, setting the real 
part of (ISBT l equal to zero leads to 



Aw x ,^(t) 



1 * 3 (£*ty) 



2 * 2 (o;^) 
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(86) 



along a ridge. Note that the leading order term in Au> x ^(t) 
is second order in St along the ridge, in agreement with 
the assumption that s lies within the instantaneous frequency 
neighborhood. In terms of scale, the ridge curve is then given 
from (gB by s(t) = [1 + Aw Xj ^(i)] oj^/uj x {t). 
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